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ABSTRACT 

We introduce a new Rigid-Field Hydrodynamics approach to modeling the magne- 
tospheres of massive stars in the limit of very-strong magnetic fields. Treating the 
field lines as effectively rigid, we develop hydrodynamical equations describing the 
1-dimensional fiow along each, subject to pressure, radiative, gravitational, and cen- 
trifugal forces. We solve these equations numerically for a large ensemble of field lines, 
to build up a 3-dimensional time-dependent simulation of a model star with param- 
eters similar to the archetypal Bp star a Ori E. Since the fiow along each field line 
can be solved for independently of other field lines, the computational cost of this 
approach is a fraction of an equivalent magnetohydrodynamical treatment. 

The simulations confirm many of the predictions of previous analytical and nu- 
merical studies. Collisions between wind streams from opposing magnetic hemispheres 
lead to strong shock heating. The post-shock plasma cools initially via X-ray emission, 
and eventually accumulates into a warped, rigidly rotating disk defined by the locus of 
minima of the effective (gravitational plus centrifugal) potential. But a number of novel 
results also emerge. For field lines extending far from the star, the rapid area diver- 
gence enhances the radiative acceleration of the wind, resulting in high shock velocities 
(up to ~ 3, OOOkms"^) and hard X-rays. Moreover, the release of centrifugal potential 
energy continues to heat the wind plasma after the shocks, up to temperatures around 
twice those achieved at the shocks themselves. Finally, in some circumstances the cool 
plasma in the accumulating disk can oscillate about its equilibrium position, possibly 
due to radiative cooling instabilities in the adjacent post-shock regions. 

Key words: hydrodynamics - stars: magnetic fields - stars: rotation - stars: mass- 
loss - X-rays: stars - gamma-rays: theory 



1 INTRODUCTION 

During their main- sequence evolution, massive, hot stars 
lack the envelope convection zones that generate magnetic 
fields in the Sun and other cool stars. Nonetheless, since 
the 1970s it has been known that a small, chemically pe- 
culiar subset - the Bp stars - possess global-scale fields at 



the kilogauss level (e.g., Borra k, Landstreet 1979 Borra, 
[Lan dstreet Thompson 1983). Moreover, the significant ad- 
vances in spectropolarimetric instrumentation over the past 
three decades have led to the discovery of ^ 100 — 1, 000 G 
fields in a number of ot her massive stars, in cluding two O- 
type stars {6^ Ori C 



2002 



HD 191612 



Donati et al 

Donati et al."2006), X-ray bright B-type stars (r Sco 
nati et al.||20Q6J, and 



Do- 



stars dHubrig et al.|2Q07 ) 

On the theoretical side 



number of slowly pulsating B-type 
the genesis of massive-star 



fields remain the subject of some controversy, with fossil- 
origin explanations (e.g., Ferrario & Wickramasinghe 2005^, 
^2006 ) competing against dynamo models involving processes 
such as core convection (Charbonneau & MacGregor 2001]), 
Tayler-Spruit instabilities ( Mullan MacDonald|20Q5( ), and 
global Rossby modes (| Airapetian||2000| ). However, consid- 
erable progress has been made in understanding how the 
magnetic fields channel and confine the stars' dense, super- 
sonic, radiatively driven winds. The seminal Magnetically 



Confined Wind Shock (MCWS) model of [Babel Mont^ 
[merle] (|1997a|b| conjectured that wind streams from op- 
posing footpoints collide near the summits of closed mag- 
netic loops, shock-heating the plasma to temperatures T ~ 
10^ - 10'^ K at which thermal X -ray emission becomes im- 
portant. Subsequent magnetohydrodynamical (MHD) sim- 
ulations by ud-Doula & Owocki (2002) confirmed the basic 
MCWS paradigm, and led to the development of a quantita- 
tive magnetic wind-shock model for the hard X-ray emission 
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of 0^ Ori C, which shows good agreement with Chandra ob- 
servations of the star ( Gagne et al.|2005p . 

MHD simulation is a powerful tool for modeling mag- 
netic wind confinement, but becomes increasingly difficult 
toward large values of the confinement parameter 77, defined 
as the ratio between magnetic and kinetic energy densities 
(ud-Doula & Owocki 2002). At large 77, the field lines are 
scarcely affected by the plasma flowing along them. Such a 
rigid character implies a high Alfven speed, and in turn a 
short numerical timestep in MHD codes to ensure Courant 



stability (Press et al. 1992). In the case of the Bp stars - 
characterized by confinement parameters up to 77 '-^ 10^ - 
the required timesteps are in fact too short for MHD simu- 
lation to be practical. 

For these strongly magnetic stars [Townsend Owocki] 



(2005 hereafter TO05) expanded on previous work by 
Michel Sturrock| (|1974) and Nakajima| ( |1985| ), to develop 
a Rigidly Rotating Magnetosphere (RRM) model based on 
the simplifying ansatz that field lines are completely rigid. 
The RRM model does not consider the detailed physics of 
the wind streams feeding into the collision shocks, but in- 
stead focuses on the fate of the post-shock plasma after it 
has radiatively cooled back down to photospheric temper- 
atures. In a rotating star, this plasma has a tendency to 
settle into magnetohydrostatic stratifications centred on lo- 
cal minima of the effective (gravitational plus centrifugal) 
potential. For an oblique dipole magnetic field, the locus of 
these potential wells resembles a warped disk that corotates 
rigidly with the star. When applied to the archetypal mag- 
netic Bp star a Ori E (HD 37479; B2Vpe), the Ha emission 
from the disk plasma shows very good agreement with that 



seen in observations ( [Townsend, Owocki Groote 2005), 
lending strong support to the model. 

Building on the success of the RRM model, this pa- 
per presents a new Rigid-Field Hydrodynamics (RFHD) 
approach to modeling massive-star magnetospheres in the 
strong- field limit. We again assume that the field lines be- 
have as completely rigid, but we now explicitly consider 
the time-dependent evolution of the magnetically channeled 
wind. This approach not only furnishes a dynamical pic- 
ture of disk accumulation, it also opens up the possibility of 
synthesizing observables for the shock-heated wind plasma, 
across a broad range of wavelengths extending from X-ray 
through to radio. 

In the following section, we consider the 1-D hydrody- 
namical problem of flow along each rigid field line, with an 
emphasis on the specific, simple case of a dipole field topol- 
ogy. In ^ we introduce a numerical code that solves the 
governing equations along many field lines, to build up a 3- 
D simulation of a massive- star magnet osphere. We use this 
code in Qto model a star loosely based on a Ori E; results 
from these simulations are presented and analyzed in ^ In 
^ we examine some of the broader issues pertaining to the 
RFHD approach, and in ^we summarize the paper. 



2 RIGID FIELD HYDRODYNAMICS 

As we discuss above, the key notion of the RFHD approach is 
that the magnetic field at sufficiently high strengths behaves 
as if it were rigid. This rigid field is anchored to the star, 
and co-rotates with it. Under the frozen flux condition of 



ideal MHD, plasma is constrained to flow along field lines, 
and it therefore describes trajectories that are fixed in the 
co-rotating frame. 

The shape of these trajectories is specified a priori by 
the chosen magnetic topology. However, the plasma state 
(density, temperature, velocity, etc.) along each field line 
is determined by the 1-D hydrodynamical problem of flow 
along a tube with changing cross-sectional area. Here, the 
'tube' can be identified explicitly with a magnetic flux tube, 
whose area varies inversely with the local magnetic flux den- 
sity B = |B| in order to ensure that V • B = 0. The char- 
acter of the flow is dictated primarily by the forces that act 
to accelerate or decelerate the plasma: pressure gradients, 
gravity, the centrifugal force, and radiative forces. Perhaps 
surprisingly, magnetic and Coriolis forces play no direct role 
in the 1-D flow problem, because they are always directed 
perpendicular to the instantaneous velocity vector v. (This 
vector is itself everywhere parallel to the field-line tangent 
vector Gs = B/B.) In fact, these forces act similarly to the 
centripetal force of a circular orbit, furnishing a net accel- 
eration perpendicular to v that leads to curved plasma tra- 
jectories yet does no work. 



2.1 Euler equations 

To elaborate on the foregoing discussion, we introduce the 
Euler equations in conservation form for the 1-D hydrody- 
namical problem comprising RFHD, 



dp Id,,. 



(1) 



^ + ^ ^ [^^ {pe + p)] = pv{gef¥ + grad) " e, + A. (3) 

Here, the independent variables are the arc distance s along 
the field line (relative to some arbitrary zero point) and time 
t, while the dependent variables are density pressure p, 
velocity = |v| and total energy per unit mass e. The term 
A{s) describes the spatially varying cross-sectional area of 
the flow tube, and depends on the magnetic topology; ex- 
pressions for this term in the case of a dipole field are de- 
rived in the following section. The acceleration vectors geff 
and grad are due to the effective gravity and the radiative 
line force, respectively, and are considered in greater detail 
in §2.3| and §2.5| Inter-relationships between p, p, v and e 
are determined from equations of state and total energy, de- 
fined in §2.6| Finally, the term A is the volumetric energy 
loss rate due to cooling processes, and is discussed in §2.7| 
For reasons elaborated there, we do not include the effects of 
thermal conduction in the energy conservation equation ([3| . 

By setting all velocities and time derivatives to zero, the 
momentum conservation equation ^ reduces to the con- 
dition of magnetohydrostatic equilibrium, in which body 
forces are balanced by pressure gradients. Because it fur- 
nishes the basis of the precursor RRM model, we review 
this static limit in Appendix [A] 



2.2 Dipole field geometry 

Although the RFHD approach is in principle applicable to 
arbitrary magnetic topologies, the present study focuses on 
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Figure 1. An illustration of the oblique dipole geometry de- 
scribed in the text. The star is shown as an oblate spheroid 
centred on the origin O, with magnetic axis M and rotational 
axis R; the angle MOR is the obliquity [3. The surfaces Sxz and 
Syz are the x-z and y-z planes of the magnetic reference frame, 
corresponding to azimuths ^ = (0°,180°) and ^ = (90°, 270°), 
respectively. Two selected field lines, having ^ = 0° and cf) = 90°, 
are shown by dotted lines; on the = 0° field line, E labels the 
magnetic equator, and P a point on the field line. The line OE 
(measured by convention in units of the rotational polar radius 
Rp) gives the magnetic shell parameter L; likewise, EP defines 
the arc distance coordinate s of the point P. (In this case, s < 0, 
because P lies in the northern magnetic hemisphere). OP is the 
radial coordinate f = r of P, and the angle MOP (ROP) gives 
the corresponding colatitude 6 (0) in the magnetic (rotational) 
reference frame. 



the simple case of an oblique dipole field. Let (r, ^,0) be the 
spherical polar coordinates in the reference frame aligned 
with the rotation axis; likewise, let (f, 0, 0) be the corre- 
sponding coordinates in the frame aligned with the dipole 
magnetic axis. (This is the same notation as adopted in 
TO05). As illustrated in Fig. [l] the magnetic axis is tilted 
with respect to the rotation axis by the magnetic obliquity 
f3. 

In the magnetic reference frame, the magnetic flux vec- 
tor is expressed as 

(^2cos6'ef + sin6'e^~) , (4) 



2(r7i?p)3 

where Bq sets the overall field strengtlj^ Rp is the rotational 

polar radius (a convenient normalizing length), and and 

are the unit basis vectors in the magnetic radial and polar 

directions, respectively. From this expression, the tangent 

vector Gs is obtained as 

B 1 / ~ ~ \ 

Gs = — = —^^=^^= 1 2 cos Oer -\- sin 6* ) . (5) 

^ V 1 + 3 cos2 e ^ ^ 

^ For an aligned dipole (f3 = 0), Bq corresponds to the polar field 
strength; however, this does not generally hold when f3 > 0, due 



This result also follows from the parametric equation for a 
dipole field line. 



r 

Rr) 



L sin 



(6) 



Nakajima|1985| [Babel k Mont merle] 1997a] ), where the 



(e.g. 

magnetic shell parameter L measures the maximal radius 
reached by the field line, in units of Rp. To label a field line 
uniquely, it suffices to specify L and the magnetic azimuthal 
coordinate (f) defining the half-plane that contains the line. 

The spatial variable s in the 1-D hydrodynamical equa- 
tions ([l]{3]) is the arc distance along each field line, and is 
found from the dipole line element 

ds^ = dr^ + dO^ = L'^rI sin (l + 3 cos^ 0^ dO'^ . (7) 

Solving this differential equation for s, we obtain 



s 

Rr> 



sinh 



+ 



cos ^ 1 + 3 cos^ ( 



(8) 



where the constant of integration is chosen to place the ori- 
gin s = at the magnetic equator, = 90°. The negative 
sign on the right-hand side arises from selecting the positive 
root of eqn. ([7|, so that s increases in the same direction 
as 0, and s is negative (positive) in the northern (southern) 
magnetic hemisphere. Since this equation is transcendental, 
calculation of the inverse function 0{s) must be undertaken 
numerically (^3.2). 



The footpoints of the field line in the northern and 
southern magnetic hemispheres are denoted as sn and ss, 
respectively. For a spherical star with radius i?*, 



R^ 



ss 
R. 



sinh-^ V^3 - 3/L 



l/Ly^I^S/l] . (9) 



However, for an oblate star (cf. ^2.4[) calculation of sn and 



Ss once again must proceed numerically (J 3.2). 

The variation in the cross-sectional area of each flow 
tube is determined by the requirement of magnetic flux con- 
servation. If Aeq is the area at the magnetic equator, then 
the area A at any colatitude must satisfy 



where 



AB ■ 



5e, 



^eq Beq , 



Bo 
2L3 



(10) 



(11) 



is the fleld strength at the magnetic equator. Eliminating B 
with the help of eqns. Q and (|6|, we flnd that 



A 

^ea 



sin^6' 



Vl + 3 cos2 



(12) 



to the oblateness of the star ({2 A) 



Combined with the inverse function 0(s), this expression is 
used to construct the area function A{s) appearing in the 
Euler equations ([l}{3]). 

In addition to the area function A, flnite-volume hydro- 
dynamical codes such as VH-1 (cf. ^ require specification 
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of the volume function 



/ 



V= Ads. 



(13) 



To obtain V in the present case, we note from eqns. ([7| 
and (IT2I that 



Ads = A^dO = LAeaRr, sin dO . (14) 
dO 



It therefore follows that 
V 



J- ("25 COS ^ - 5 cos W + cos 5^) + 
320 V J 



448 



cos 70 



, (15) 



and, as with A{s), the inverse function 0{s) is used to find 
V{s). 

2.3 Effective potential 

The effective gravitational acceleration geff in eqns. ( 2|3 ) 
combines the Newtonian gravity with the centrifugal force 
associated with enforced co-rotation. Together, these forces 
are derived from a scalar effective potential $eff, 

geff = -V$eff . (16) 

Within the Roche (point-mass) approximation, this effective 
potential is given by 



. GM. 1 2 2 

^eff = -il TaJ 

r 2 



(17) 



Here, M* is the stellar mass, Q the angular rotation fre- 
quency, and uj = rsin^ is the distance from the rotation 
axis. Combining the above two expressions, the effective 
gravity term in the momentum and energy conservation 
equations (2[3) is found as 

geff • es = ■ es + il tue^ • es , (18) 

where Gr and e^ are the unit vectors in the r and ui direc- 



tions, respectively. For a dipole field (^2.2), 



ef • e^ 



2 cos 6* 



^ _ cos X (19) 

VI + 3cos2 6' 

gives the component of the local radial vector projected 
along the field line. The angle x introduced in this expres- 
sion, being that between the field line and the radial vector, 
also appears below in the equations governing the radiative 
acceleration. 

2.4 Stellar properties 

In addition to generating the effective gravity, $eff de- 
termines the shape and surface properties of the oblate, 
centrifugally distorted star. In the Roche approximation 



(eqn. 17), the surface is an equipotential whose radius r* 



varies with rotational colatitude as 



Rry 



w sin 

(e.g., Cranmer]|1996 ). Here, 



TT + COS ^ (w sin 0) 



8GM, 



(20) 



(21) 



is the normalized rotation angular frequency, with w — 1 
corresponding to critical rotation. 

Due to gravity darkening (von Zeipel 1924), the flux 



emitted by a rotating, radiative stellar envelope varies in 
proportion to the local effective gravity |geff|. However, in 



evaluating the radiative acceleration (^2.5) and the rate of 
inverse Compton cooling (J 2.6), we choose for simplicity to 
treat the circumstellar radiation field as originating from 
a spherically-symmetric point source of luminosity L* . For 
consistency, we therefore neglect the variation of the surface 
flux, setting it to the constant value 



(22) 



Here, So is the total surface area of the oblate star, which 
can be approximated to within 2% by the polynomial 

So = 47ri?p (1 + 0.19444^^ + 0.28053^^ - 1.9014^^ + 

6.8298^^ - 9.5002^/;^° + 4.6631^^^) (23) 



( |Cranmer|1996 ). From and the Stefan-Boltzmann law, a 
nominal stellar surface temperature is defined as 



2.5 Radiative acceleration 



1/4 



(24) 



To obtain an expression for the radiative acceleration grad 
we employ the Castor, Abbott & Klein (1975 hereafter 



CAK) formalism for line-driven stellar winds. For simplicity, 
the star is treated as a point source of radiation, giving an 
acceleration 



grad 



1 KeL^Q / l^^l 



1 — a 47rr^c \pcQK,e 



(25) 



(e.g., Owocki 2004|). Here, a is the CAK power-law index; 
Q is the dimensionless line strength parameter introduced 
is the stellar luminosity (^2.4); and Ke 

for a fully 



by.Gaylev (1995) 



is the electron-scattering opacity, ^ 0.34 cm 
ionized solar-composition plasma. The term 6 V , a measure 
of the local velocity gradient that determines the |Sobolev| 
( 1960| optical depth, is defined as 



From eqn. ( 19 ) 



and thus 



6v = Gr ■ V(er • v) . 



V COS X ? 



Sv = er ■ V{v COS x) — COS X • . 



(26) 



(27) 



(28) 



(The second equality follows because x is a function of 
alone, and therefore commutes with the radial directional 
derivative er • V). Expanding out the gradient operator in 
the magnetic reference frame, we obtain 



6v — cos X 



(29) 



To evaluate the radial derivative in this latter equation, 
it is necessary to know the flow velocity on field lines adja- 
cent to the one under consideration. In order to avoid this 
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complication, we make the simplifying assumption that the 
polar velocity derivate vanishes, 



! 0. 



(30) 



As we demonstrate in Appendix [B] this assumption implies 
that 



dv\ f dv 



dr I Q \dsy^ 



secx • 



Substituting this expression back into eqn. ( 29 ) yields 



(31) 



(32) 



which in combination with eqn. (25) leads to the final ex- 
pression for the radiative acceleration, 

_ 1 KeL.Q f \dv/ds\ 
^^^^ 1 — a Anr'^c \ pcQKe 



(33) 



To conform with our earlier notation (cf. ^ 2.1 ), here we have 
dropped the 'L' subscript on the velocity gradient dv/ds. 

The central approximation here, eqn. (30), is equivalent 
to assuming that v - which might more correctly be termed 
the flow speed - depends only on f, although of course the 
local flow direction still depends on since it follows the 
field-line orientation. Analysis based on 2-D MHD simula- 



proximation is well-justified in the regions of rapid wind ac- 
celeration near to the star. Moreover, in the present work 
we have carried out a posteriori checks using an ensemble of 
neighboring 1-D RFHD calculations, and find that the error 
in grad arising from using eqn. (32) for Sv instead of the ex- 



act expression (29) is typically 10% or less throughout the 
radiatively driven regions of the magnetosphere. 

We emphasize that the benefit from this approximation 



is significant. Because eqn. (33) depends only on the velocity 
gradient along a field line, it allows the flow to be modeled 
completely independently of other field lines. Not only is 
this advantageous in terms of computational efficiency, it 
also simplifies the interpretation and analysis of simulation 
results (^. 



2.6 Equations of state and energy 

We assume an ideal gas, so that 

pkT 

P , 

pu 



(34) 



with T the temperature and u the atomic mass unit. The 
mean molecular weight fx is determined from an expression 
appropriate to a fully ionized mixture. 



2X+-(1. 



■X-Z) + 



(35) 



where X and Z are the usual hydrogen and metal mass 
fractions. The accompanying equation for the specific (per- 
unit-mass) total energy is 



1 kT 
7 — 1 pu 2 

where 7 is the ratio of specific heats. 



(36) 



2.7 Cooling and thermal conduction 

The volumetric cooling rate A in the energy conservation 
equation ^ is evaluated as the sum of two terms. 



A = Aat + Aic , 



(37) 



representing contributions from atomic processes and in- 
verse Compton scattering by thermal electrons, respectively. 
The first term is calculated from 



where 



and 



Aat = nenpC{T) , 



X 



1±X_ 
2u ^ 



(38) 



(39) 



(40) 



define the proton and electron number densities, respec- 
tively. The temperature-dependent cooling coefficient C is 
obtained from the curve published by |MacDonald Bailey] 
((1981). 

The inverse Compton term Aic is evaluated with the aid 



of eqn. (4) of White k Chen (1995), 



Aic 



4—neUphkT. 
c 



(41) 



(Note that in their eqn. 5, these authors define the symbol 
Aic differently than above). Here, 



f7ph 



(42) 



is the energy density associated with the star's radiation 
field, evaluated in the same point-star limit as the radiative 



acceleration (J 2.5). 

We have elected to neglect the effects of thermal conduc- 
tion in the energy conservation equation, because to include 
these properly requires addressing a large number of issues 
that are beyond the present scope. These center on uncer- 
tainties both in the treatment of departures from the clas- 



sical Spitzer 


( 1962 ) thermal conductivity (e. 


Eichler|1992 


Pistinner k Eichler|1998j), and 



ration of heat flux saturation near shocks and other regions 
of steep temperature gradients (see Lacey 1988| Bandiera & 



Chen |1994 , and references therein). Moreover, apart from 



these physical difficulties, proper inclusion of conduction in 
a hydrodynamical code is a challenging task, particularly in 
the context of accurately modeling the kind of strong shocks 
that occur in our simulations (see, e.g., |Reale|1995t . We thus 



defer consideration of thermal conduction to future work. 



3 THE RFHD CODE 



As discussed in §2.5| the utility of the approximate expres- 
sion (32) for the velocity gradient term Sv lies in the fact 



that the fiow along each individual field line may be sim- 
ulated completely independently of other field lines. Thus, 
by performing many separate 1-D simulations for differing 
field lines and piecing them together, a 3-D hydrodynami- 
cal metasimulation of a massive-star magnetosphere can be 
built up at a fraction of the computational cost of an equiv- 
alent MHD calculation. 
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Each individual 1-D simulation requires solution of the 
Euler equations ([l]{3|, and for this we employ a customized 
version of the VH-1 hydrodynamical code developed by J. 
Blondin and colleagues. VH-1 is a finite volume code based 
on the Lagrangian version of the piecewise parabolic method 
(PPM) devised by |Colella &: Woodward] ( |1984| ). The modi- 
fications to VH-1 primarily encompass incorporation of the 
dipole geometry (J 2.2), acceleration terms (§^ 2.3|2.5 ) and 
cooling (^2.7) specific to the RFHD problem. We do not 
discuss these modifications in detail, but in the following 
sections we highlight specific issues that arose during the 
code development phase. 



3.1 Grid design 

The Eulerian grid in VH-1 must be designed with care, to en- 
sure that regions of physical interest are properly resolved, 
and to avoid the generation of spurious numerical instabili- 
ties. To this end, we divide the grid into three domains. The 
two 'surface' domains are each composed of a fixed number 
A^s of zones, extending from the surface footpoints s = sn, ss 



(cf. ^3.2) to somewhat above the expected position of the 
sonic point \v\ = a*. (Here a* = y/kJ^JJm denotes the 
isothermal sound speed at the surface, with T* the stellar 
surface temperature introduced in p.4[ ) These zones are 
non- uniform, with the size of each zone being 1.1 times that 
of its neighbour closer to the stellar surface. With their high 
spatial resolution, the surface domains are designed to re- 
solve the smooth transition of the wind from a subsonic, 
near-hydrostatic state to a supersonic outflow. 

The third, ' magnet osphere' domain spans the physically 
interesting regions of the circumstellar environment, and is 
composed of A^m zones of uniform size ds ^ (ss — sn)/A^s, 
that bridge between the two thin surface domains. In deter- 
mining an appropriate value for A^m, we are motivated by a 
desire to resolve the dense, co-rotating disk predicted by the 
RRM model to accumulate at local minima of the effective 
potential (see TO05; see also Appendix[A|. As demonstrated 
in the Appendix, the scale height of this disk is independent 
of the magnetic shell parameter L when L becomes large. 
Thus, we vary A^m with L according to the simple formula 



int(Adiv + A*L' 



-l/2x 



(43) 



where int() denotes the integer part. The first term in the 
parentheses ensures that the disk remains properly resolved 
at large L, with a constant number of zones per asymptotic 
scale height hoo (cf. eqn. A9); while the second is an ad- 



hoc addition that increases the spatial resolution near the 
star. The parameters Ad and A* are adjusted to achieve the 
desired spatial resolution in the far- and near-star limits. 



3.2 Field-line coordinates 

To determine the arc distance coordinates (sn,ss) of the 
field-line footpoints, required in setting up the Eulerian grid 



(^3.1 ), we solve the equation 



LRp sin 



(44) 



to find the magnetic colatitudes (^n,^s) associated with 
each footpoint. (The footpoints are not mirror symmetric 
about the magnetic equator unless i(; = or = 0°). The 



corresponding values of s then follow from eqn. (|8|. In the 
above equation, is a function of the rotational colatitude 
(cf. eqn. 20 ); thus, it depends implicitly on and via the 
coordinate transformation between magnetic and rotational 
reference frames (s ee Fig.p]). To solv e this equation, we use 
Brent's algorithm ( [Press et al.|1992| ). 

Once the Eulerian grid is established, we calculate 
at each zone boundary from the corresponding s values, by 
inverting eqn. ([s]). For this, we use a Newton- Raphson iter- 
ation operating on cos^, starting from the initial guess 



LRry 



sinh ^ 
2a/3 



+ 1 



(45) 



After each Lagrangian step of the PPM algorithm, the ad- 
vection of zone boundaries means that the values must be 
updated. The Newton- Raphson iteration is in this case too 
slow to be useful, so we instead evaluate by cubic-spline in- 
terpolation ( Press et al.|1992 ) from the stored Eulerian-grid 
values. 



3.3 Initial state 

The initial, t = state for simulations is based on a spheri- 
cally symmetric, accelerating wind. The velocity is obtained 
by projecting a canonical velocity law onto field lines: 



V = Voo ^ - cosx • 



(46) 



Here, Voo is the terminal velocity obtained from a non- 
rotating, point-star CAK model. 



2GM^ 



1 — a Rp 



(47) 



The density is likewise obtained from the condition of steady 
spherical outflow. 



where now 



M = 



a 



M 



QFe 



a — 1 

is the CAK mass-loss rate, with 



(l-a)/a 



(48) 



(49) 



(50) 



the Eddington parameter associated with electron scatter- 
ing. By assuming an isothermal initial flow, with T = T*, 
the pressure p is obtained from the equation of state ( 34 ) . 



3.4 Boundary conditions 

Boundary conditions are implemented in VH-1 by setting 
dependent variables (p, p, e) in ghost zones at both ends 
of the grid. To allow for the possibility of outflow and inflow, 
we fix only the density and pressure in the ghost zones, as 
p = and p = p^. The stellar density can be chosen with 
some degree of latitude, so long as it remains appreciably 
above the wind density 



M 



(51) 
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(Mq) 



Rp 



(K) 



(L0) 



(3 

n 



X 



5.3 



22,500 7,420 55 0.7 0.02 



^rot 

(days) 



7 



A* 



1.2 500 0.6 5/3 25 100 100 
Table 1. Parameters for the RFHD simulation described in ^ 

at the sonic point; we find that a choice = 10 ps gives a 
smooth and steady wind outflow. The corresponding stellar 
pressure is obtained from the equation of state ( 34 ) , with 
and T 



p — and 1 = T*. The velocity in the ghost zones is 
linearly extrapolated from the first pair of computational 
zones (subject to the constraint that \v\ does not exceed 
a*), and the total energy per unit mass is evaluated using 
eqn. (36). 



3.5 Cooling 

Cooling in the customized VH-1 is time-split from the rest 
of the PPM algorithm. Across each timestep dt, a new tem- 
perature Tj in the j'th zone is calculated from the current 
temperature Tj and density pj via 



T] = T, + 



{^-l)f,uKT},Pj) + KTj,Pj) 



k 



2pj 



(52) 



where A(T, p) is the cooling rate discussed in J 2.7 Since this 
equation is implicit, the code uses 3-step iteration to con- 
verge toward an approximate value for T|. In zones where 
the temperature would drop below the stellar surface tem- 
perature T*, it is reset to T*; this reflects the tendency for 
the photospheric radiation field to keep circumstellar plasma 
warm. The resulting T| values are used to update the pres- 



sure in each zone, through the equation of state (34). 



lines, arranged on a grid of 60 magnetic shell parameters 
L and 19 magnetic azimuths <j). The azimuth values span 
the 0° ^ ^ 90° quadrant in 5° increments; symmetry is 
used to replicate the simulation results over the remaining 
90° ^ ^ 360° interval, effectively leading to 72 azimuthal 
points. The 60 values of the shell parameter are distributed 
according to the formula 



InL = lnl.2 + ln(11.2/1.2)- 



59 



(^ = 1 . . . 60) . (53) 



This logarithmic distribution provides higher resolution 
close to the star, where we expect magnetospheric struc- 
ture to depend more sensitively on L. The innermost, I — 1 
field lines have L = 1.2; by comparison, the stellar equatorial 
radius for the adopted parameters is i^eq = 1.11 i?p. These 
field lines are composed of A^s = 25 grid zones in each of 
the two surface domains and - with the choices Ad and A* 
listed in Table [l]- AT^ = 211 zones in the intervening mag- 
netosphere domain. The outermost, ^ = 60 field lines have 
L — 11.2, and therefore extend out to just over six times 
the Kepler co-rotation radius tk = 1.86 i?p (cf. eqn. A8). 
These field lines again have A^s = 25, but now the magne- 
tosphere domain has A^m = 1149 zones, ensuring a zone size 
ds = 0.027 i?p that is less than half the asymptotic scale 
height hoo = 0.064 i?p defined by eqn. (|A9|. 



Each 1-D simulation runs for 200 stellar rotation cy- 
cles (~ 21 Ms). This is significantly longer than typical flow 
times 20 ks, to ensure that radiative relaxation from the 
initial state is obtained, and to follow the long-term evo- 
lution of the magnetosphere. Since the simulations for each 
field line are independent, they can be distributed across dif- 
ferent processors, computers or even clusters. In the present 
case, the 1140 individual simulations were run on a cluster 
of eight 2.0 GHz dual-core AMD Opteron nodes (for a total 
of 16 processors). On average, each simulation took ~ 200 
minutes (although there was a significant spread about this 
mean), giving a total computation time of ^ 10 days. 



4 CALCULATIONS 

As an initial test of the code discussed in the preceding sec- 
tions, and to explore the RFHD approach in general, we 
develop a model for the magnetosphere of an oblique-dipole 
star whose parameters (Table [T]) loosely coincide with those 
of cr Ori E. (We stress that we are not attempting to fine- 
tune a model for a Ori E; however, it makes sense to begin 
our qualitative investigations in the same approximate re- 
gion of parameter space as this archetypal star.) The mass, 
radius and surface temperature of the star are taken from 



the recent study by Krticka et al. (2006), with the luminos- 



ity calculated from eqns. (22) and (24). The wind parame- 



ters Q and a are assigned values typical to hot stars (e.g., 
Gayley|[l995 ), resulting in a CAK mass loss rate (eqn. 49) 



of 3.7 X 10 ^ Mq yr ^. The rotation period Prot is based 



on the photometric period measured by |Hesser, Ugarte &| 
Moreno ( 1977), and leads to a normalized rotation frequency 



w — 0.73. The dipole obliquity {3 is adopted from the RRM 



model presented by Townsend et al. (2005). 

To piece together a 3-D model for the magnetosphere, 
we use VH-1 to simulate the flow along a total of 1140 field 



5 RESULTS 

The complete 3-D magnetospheric model we describe above, 
composed of 1140 1-D simulations, occupies ^ 36 Gb of stor- 
age. Presenting this significant dataset in an informative 
manner poses quite a challenge. Fortunately, the indepen- 
dence of the individual simulations is once again of benefit, 
since the flow can be analyzed first in one dimension, along 



an individual field line (^5.1 ); then in two dimensions, along 
field lines lying in the same meridional plane (J 5.2); and 



then in three dimensions (^5.3), for the complete model. 



5.1 1-D, along individual field line 

In this section, we focus exclusively on the flow along the 
field line having a magnetic shell index ^ = 20 (cf. eqn. 53), 
corresponding to L = 2.45, and an azimuth = 90°. With 
these parameters, the field line extends out beyond the Ke- 
pler radius tk = 1.86 i?p, and passes through the intersec- 
tion between magnetic and rotational equators. According 
to the RRM model (cf. TO05), these properties should be 
favourable to the steady accumulation of cooled plasma at 
the field- line summit. This expectation is amply confirmed 
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log p (g cm"^) log T (K) V (10^ km s"^) 



-16 -15 -14 -13 -12 45678 -2 -1012 




-2-1012 -2-1012 -2-1012 

5 (^p) 5 (^p) 5 (^p) 



Figure 2. The time evolution of the flow along the (L, ^) = (2.45, 90°) fleld line, showing the density p, temperature T, and velocity i; as a 
function of arc distance s and time t. The magnetic equator is situated at s = 0, and the two footpoints at s = (sn, 5s) = (—2.28, 2.28) R^. 
The horizontal dotted lines, labeled alphabetically, show the locations of the snapshots plotted in Fig.js] the 'A' line is situated at t = Os. 




-2-1012 -2-1012 -2-1012 

5 (^p) s (^p) 5 (^p) 

Figure 3. Snapshots of the flow along the (L, ^) = (2.45,90°) fleld line, plotting the density p, temperature T, and velocity as a 
function of arc distance s. The times of the snapshots are indicated in Fig. [2] 
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by Fig. [2] which shows the time evolution of p, T, and v 
along the field line for the first 3 Ms of the simulation. (The 
dynamics during the remainder of the simulation are not 
significantly different from those seen toward the end of this 
initial timespan) . The figure should be interpreted with the 
aid of Fig. |3] which plots snapshots of the flow variables at 
four epochs of interest. 

The first, 'A' snapshot shows the initial configuration 
described in §3.3| Because the velocity in both hemispheres 
is supersonic, reverse shocks quickly form at the magnetic 
equator s = 0, and then proceed to propagate back down 
the field line toward each footpoint. As wind plasma flows 
through these shocks it experiences an abrupt reduction in 
velocity, matched by corresponding discontinuous increases 
in density and temperature. The temperature jump, from 
T = T* to T ^ 2 X 10^ K, mean that the post-shock plasma 
cools initially with photon energies kT ^ 2 keV in the X-ray 
range. 

The downward propagation of the shocks halts when the 
ram pressure of the pre-shock plasma matches the gas pres- 
sure in the post-shock regions. The resulting quasi-steady 
state, composed of standing shocks at s = ±1.6 i^p that en- 
close hot, post-shock cooling regions, is shown in snapshot 
'B'. Situated at the centre of these regions is the cooled 
plasma predicted by the RRM model. This plasma accu- 
mulates at a local minimum of the effective potential $eff 
(as sampled along field lines), where it is supported in sta- 
ble magnetohydrostatic equilibrium by the centrifugal force 
(see TO05; see also Appendix [A|. As we discuss further in 
§5.3| the locus formed by such potential minima resembles 
an azimuthally warped disk. 

To characterize the behaviour of the cool disk plasma, 
we briefly digress to introduce the three moments 



CTd 



Sd = 



J psA ds , 



and 



hd 



^eqCTd 



^eqCTd 



I Pis 



Sd) Ads 



1/2 



(54) 
(55) 

(56) 



where the integrals extend over the disk regions having T = 
T*. These moments represent, respectively, the disk surface 
density (projected into the magnetic equatorial plane), the 
disk centroid, and the disk scale height. When applied to 
the 1-D simulations, the above expressions are evaluated via 
finite- volume equivalents. 



CTd = ^ ^ Pj dVj , 



hd = 



-4eqCrd 



SdYdVi 



1/2 



(57) 
(58) 

(59) 



where dV^ is the volume of the j'th zone, and sj the arc 
distance of the zone's centre. As before, the summations 
extend over disk zones having Tj = . 

Fig. [i] plots the moments (57 _59) for the (L, 0) = 
(2.45, 90°) field line. The rightmost panel indicates that dur- 
ing the early stages of the simulation, the disk thickness is 



substantially smaller than the value ho = 0.09 i?p predicted 
by the RRM model (cf. eqn.[A6|. Indeed, up until t = 0.6 Ms 



hd is identically zero, indicating that the disk extends over 
only a single zone. The reason for this confinement is that 
the internal pressure of the cool disk is insufficient to support 
it against the relatively high pressure of the hot plasma in 
the surrounding post-shock regions; thus, the disk becomes 
compressed into the smallest volume resolvable by VH-1. 

This situation changes once the steady accumulation of 
plasma (Fig. |4] leftmost panel) raises the internal pressure 
to a point where the disk suddenly grows to encompass a 
greater number of zones. Three of these expansion events 
are apparent in Fig. [2] at t = 0.6, 1.6 and 2.1 Ms. After each 
event, the disk plasma undergoes oscillations back and forth 
about its equilibrium position s = 0; these are most readily 
seen in the centre panel of Fig. ^ for t > 2.1 Ms. The os- 
cillations are accompanied by variability of the post-shock 
cooling regions, involving the abrupt movement of the shock 
fronts toward s = 0, followed by a rebuilding of the cool- 
ing regions via wind feeding. Snapshot 'C in Fig. |3] shows 
the state of the flow immediately after an ingress of the 
northern-hemisphere (s < 0) shock; observe how this shock 
is 0.45 i?p closer to the disk than in snapshot 'B'. Because 
the density jump across a strong shock scales proportionally 
to the pre-shock density, this movement of the shock front 
translates into a reduction in the post-shock density. Thus, 
p ^ 1.3 X 10~^^gcm~^ just downstream of the northern 
shock in snapshot 'B', but p ~ 0.4 x lO^^^gcm""^ down- 
stream of the shock in snapshot 'C. 

The disk oscillations after the first expansion event are 
damped, dying out over a timescale ^0.5 Ms. The damping 
is even stronger for the very weak oscillations seen after the 
second event. Following the third event at t = 2.1 Ms, how- 
ever, the oscillations persist at a relatively high amplitude 
all the way to the end of the simulation. Snapshot 'D' in 
Fig. [3] illustrates the flow when the oscillating disk plasma 
is undergoing a southward displacement (sd = 0.02 i?p). 

To explore the nature of the oscillations. Appendix |C] 
develops a linear analysis of the response of the disk plasma 
to small-amplitude departures from stationary equilibrium. 
For perturbations that conserve dd, the analysis reveals a 
spectrum of normal modes having periods 



LOrn 



^ho 
m a* 



(m=l,2,3. 



(60) 



The oscillations seen in Figs. |2j|4] are wholly consistent with 
the excitation of the dipole (m = 1) mode. In particular, the 
period Pi = 85.5 ks predicted by the expression above is in 
good agreement with the value P = 79.2 ks measured from 
the strongest peak in the Fourier transform of the Sd data 
for t > 2.1 Ms. 

What excites the dipole oscillations? In the simulations, 
the coincidence between the expansion events and the on- 
set of oscillations reveals that these events impart an ini- 
tial 'kick' to the disk plasma. (The kick originates because 
the perturbation introduced by an expansion event is usu- 
ally asymmetric.) However, some other process is clearly re- 
sponsible for maintaining and even amplifying the oscilla- 
tions, as seen after the third expansion event at t = 2.1 Ms. 
A likely candidate for this process is the cooling instabil- 
ity discussed by Chevalier & Imamura (1982). These au- 



thors show that radiative cooling of thermal plasma in the 
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Figure 4. The surface density (Jd, centroid and scale height of the disk for the (L, ^) = (2.45, 90°) field line, plotted as a function 
of time t. The dotted lines show the corresponding predictions of the RRM model. 
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Figure 5. The density p in the vicinity of the cool disk, plotted 
as a function of arc coordinate s for the (L, cj)) = (2.45, 90°) field 
line. The solid line indicates the result from the RFHD simulation 
at a time t ~ 20.2 Ms, while the dotted curve shows the prediction 
of the RRM model at the same t. 



temperature range ^ 10^ — 10^ K is linearly overstable due 
to the negative temperature exponent of the cooling coef- 
ficient C. For the parameters of the cooling regions of the 
(L,0) = (2.45,90°) field line, the periods of the fundamen- 
tal and first-overtone unstable cooling modes are P ~ 160 ks 
and P 40ks respectivel}j^ bracketing the disk oscillation 
period P = 79.2 ks. This lends support to the hypothesis 
that a coupling between the disk modes and the cooling 
modes results in a global overstability, driving the disk os- 
cillations until non-linearity leads to saturation. We intend 
to investigate this hypothesis further in a future paper. 

To bring the present section to a close, we compare 
the RFHD simulation for the (L, 0) = (2.45,90°) field line 
against the predictions of the RRM model. For the parame- 
ters specified in Table [l] the RRM modej^ predicts a rate of 



change dd = 1.38 x 10~^ gcm~^ s~^ for the disk surface den- 
sity. Given the approximations employed in developing the 
model, this value is in remarkably good agreement with the 



empirical value fid = 1.31 x 10 



"^gcm ^ 



s derived from a 



linear least-squares fit to the simulation data. 

To compare the distribution of disk plasma, Fig.[5]plots 
the density as a function of s for both model and simulation, 
at a time t ~ 20.2 Ms near the end of the simulation cho- 
sen so that the disk displacement Sd is close to zero. Once 
again, there is encouraging agreement between the RRM 
prediction and the simulation result. The only significant 
difference is that the wings of the Gaussian density distri- 



bution (cf. eqn. A5 ) are truncated in the simulation, due to 



the pressure of the hot plasma in the adjacent post-shock 
regions. 



5.2 2-D, in meridional planes 

We now extend the analysis to two dimensions, by consider- 
ing 1-D RFHD simulations all lying in the same meridional 
plane. Fig. [6] shows 2-D images of the state of the flow at 
the end of the simulations, for the cj) — 90° meridional plane. 
(The figure also shows images for the = 0° plane, but we 
defer discussion of these data until later). The cool disk of 
accumulated plasma is clearly seen along the y axis in the 
= 90° plane, surrounded by hot (T - lO'^ - 10^ K) post- 
shock cooling regions. The disk does not extend all the way 
to the star, but is instead truncated at a radius r = 1.8 i?p. 
Inside this radius, the centrifugal force is not strong enough 
to support plasma against the inward pull of gravity, and no 
accumulation occurs. Although dense knots of plasma are 
formed close to the star by compression between opposing 
wind streams, these knots quickly slide down the magnetic 
field toward one or the other of the footpoints. 

Figure [7| illustrates this fallback process in a sequence of 
density snapshots, showing the evolution of the inner parts 



^ These val ues are obtained from the a = —1 eigenfrequencies 
tabulated by [Chevalier Imamura| ( |l982| , with nin = 900 km s"-*^ 
and XsQ = —l.QRp. 

^ See TO05, their eqn. (34); in evaluating the stellar-surface field 



tilt /i* appearing in this equation (and in other expressions from 
the RRM formalism), we take into account the stellar oblateness 
due to rotation, which for simplicity was neglected in the TO05 
treatment. 
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Figure 6. The state of the flow at the end of the simulations, showing the density p, temperature T, and velocity v in the cf) = 90° 
(upper row) and ^ = 0° (lower row) meridional planes. (See Fig. ^ for a recapitulation of the geometry of these planes). The oblate 
star is indicated in grey at the left-hand side of each image. Outside the regions threaded by simulation fleld lines (i.e., for L < 1.2 and 
L > 11.2), the images are left blank. The dotted rectangle in the upper left image indicates the region shown in detail in Fig.[7| 



of the magnetosphere shortly before the end of the simula- 
tions. The white arrow in the leftmost panel indicates a knot 



that has formed at {y, 



(1.6, 0.1) Rp. In the centre panel, 



the knot has migrated further into the northern magnetic 
hemisphere, and in the rightmost panel begun its descent to 
the stellar surface. This figure is reminiscent of the infalling 
plasma seen in the MHD simulations by |ud-Doula Owockij 
(^002, their fig. 4); in fact, the only phenomenological dif- 
ference lies in the scale of the knots. Cross- field coupling in 
the MHD simulations allows coherence between the flow on 
adjacent field lines. However, this coupling is absent in the 
RFHD simulations (due to the neglect of the polar velocity 
derivative when calculating the radiative acceleration; see 
^2.5), and the scale of the knots is set therefore by the L 



grid spacing (cf. eqn. 53). 

The absence of cross-field coupling is also ultimately re- 
sponsible for the significant amount of structure seen in the 
post-shock cooling regions. If disk oscillations are excited 
on one particular field line, but are absent from an adjacent 
field line, then discontinuities in the cross- field direction ap- 
pear in the flow properties, resembling a reversed letter 'C. 
Thus, for instance, the (L, 0) = (3.88, 90°) field line in Fig.ji] 
has a significantly lower density than its neighbours because 
disk oscillations are excited on this field line, but not on the 



adjacent ones. As with the fallback, we expect in reality that 
cross- field coupling will tend to smear out such sharp discon- 
tinuities. One possible approach to including this coupling 
in RFHD simulations is discussed in §6.3. 1| 

Turning now to the temperature data in Fig. [6] a gradi- 
ent can be seen across the post-shock regions, with the outer 
parts being up to an order of magnitude hotter (T ^ 10^ K) 
than those near the star (T - lO'^K). This gradient arises 
from two distinct effects. First, field lines having larger mag- 
netic shell parameter L undergo a greater area divergence, 
and are therefore characterized by lower plasma densities 
and faster flow velocities (see |Owocki ud-Doula|2004| for a 
more in-depth discussion of this effect) . This naturally leads 
to a tendency for larger-L field lines to experience hotter 
post-shock temperatures. 

Second, in a process first conjectured by |Babel &: Mont^ 



merle ( 1997a), the post-shock plasma can be pushed to even 



higher temperatures by the action of the centrifugal force. 
To accumulate onto the cool disk, the plasma must first de- 



scend to the bottom of the effective potential well (^5.1). 
The consequent release of centrifugal potential energy heats 
the plasma, by an amount that ultimately depends on the 
distance from the rotation axis. In the (j) — 90° data shown in 
Fig.|6]this effect raises the plasma temperature on the outer, 
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Figure 7. Snapshots of the density p in the inner parts of the ^ = 90° meridional plane (see the dotted rectangle in Fig.jeJ, at a time 
t ^ 20.4 Ms near the end of the simulations (left), and then at increments of one quarter (centre) and one half (right) of a rotation cycle 
later. The white arrows indicate the location of the dense knot discussed in the text. As in Fig. [6] the star is indicated in grey at the 
left-hand side of each image. 



L = 11.2 field line from T = 1.1 x 10^ K at the shocks, to 
T = 1.8 X 10^ K adjacent to the disk. 

The centrifugal heating in the post-shock regions com- 
petes against cooling by atomic and inverse Compton pro- 
cesses (cf. J 3.5). Initially, the centrifugal effect dominates. 



because plasma densities are so low that atomic cooling (a 
process; see eqn. 38) is exceedingly inefficient. However, this 
situation is subsequently reversed; as plasma moves down- 
stream of the shocks, the density eventually becomes high 
enough for atomic processes to cool it rapidly, in thin lay- 
ers on either side of the disk (see ^5.3). During this pro- 



cess, inverse Compton cooling is relatively unimportant; at 
low densities it is more efficient than atomic cooling, but it 
never becomes the dominant term on the right-hand side of 
the energy conservation equation 

Although the foregoing discussion focuses on the flow 
in the = 90° meridional plane, it generally applies to 
other meridional planes. However, one significant exception 
concerns the location of the cool disk. The = 0° images 
in Fig. [6] reveal that the disk is not symmetric about the 
magnetic axis, as one might presume from considering the 
(f) — 90° images on their own. Rather, the disk is warped in 
the azimuthal direction. To explain this notion, let Od de- 
note the centroid magnetic colatitude of the disk, obtained 
by setting s = Sd in eqn. ([s]) and solving for 6. Then, in a 
given meridional plane the disk lies approximately along a 
straight ray emanating from the origin, with 6d remaining 
constant as L is varied. However, Od changes from plane to 
plane, resulting in azimuthal warping; for instance, Od ~ 70° 
in the = 0° meridional plane, whereas Od ~ 90° for the 
(j) = 90° case. This warping is one of the key predictions of 
the RRM model, and its origin is discussed in greater detail 
by TO05. 

In addition to revealing the disk warping, the = 0° 
images in Fig. [6] also illustrate a novel flow phenomenon. On 
a number of the innermost (L < 1.8) field lines, the northern 
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Figure 8. The disk surface density era in the = 90° meridional 
plane, plotted as a function of magnetic shell parameter L. The 
solid line shows the final state at the end of the RFHD simula- 
tions, while the dotted line indicates the corresponding prediction 
of the RRM model. 



reverse shock propagates back down to the stellar surface. 
This leads to a siphon configuration, in which plasma flows 
unidirectionally from the southern magnetic footpoint to the 
corresponding northern footpoint. (The direction of the flow 
is set by the relative location of the rotation axis; in the = 
180° meridional plane, north-to-south siphon flows occur). 
Because the southern reverse shock remains intact, part of 
the siphon flow is composed of shock-heated plasma at T ^ 
5 X 10^ K; this plasma is relatively dense due to its proximity 
to the star, and therefore makes a significant contribution to 



the soft X-ray emission from the magnetosphere (see J 5.3). 

As with the 1-D analysis in §5.1| we bring the present 
section to a close by comparing the RFHD simulation results 
against the predictions of the corresponding RRM model. 
Fig.jsjplots the disk surface density ad in the = 90° merid- 
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ional plane, at the end of the simulations and for the model. 
The agreement between the two is once again very encour- 
aging, especially in the innermost regions of the disk. In the 
outer regions, the simulations predict a rather smaller sur- 
face density (- 10-25%) than the RRM model. This modest 
discrepancy appears correlated with the observation that, 
for field lines having L > 6, the disk remains confined to a 
single zone - no expansion events occur on these field lines 
over the duration of the simulations. 



5.3 Full 3-D model 

Having examined the results of the RFHD simulations for 
an individual field line, and for field lines lying in meridional 
planes, we now turn to the complete, three-dimensional pic- 
ture. Fig. [9] shows images of the proton column density A^p 
at the end of the simulations, viewed from six equally spaced 
rotational azimuths spanning the range 0° ^ ^ 180° (the 
remaining interval 180° ^ ^ 360° is mirror-symmetric 
through the vertical axis). Following [Townsend et al. (2005 ), 



a viewing inclination i = 75° with respect to the rotation 
axis is adopted. The column density is calculated from ray 
integrals of the form 



dzi . 



(61) 



where Zi is the perpendicular distance to the image plane. 



and the proton number density Up was defined in eqn. (39). 
The mass density p appearing in this latter equation is cal- 
culated using trilinear interpolation from the [s^L^cj)) field- 
line coordinate system onto the (xi, 2/1,2:1) Cartesian image 
coordinate system. Along rays that intersect the star, the 
integral is truncated at the stellar surface. 

The figure clearly illustrates the three-dimensional 
structure of the dense, co-rotating disk discussed in previous 
sections. (The plasma in the wind and post-shock regions is 
effectively invisible, since it is orders-of-magnitude less dense 
than that comprising the disk.) This disk possesses three im- 
portant characteristics predicted by the RRM model. First, 
it has an average inclination that lies somewhere between 
the magnetic and rotational equatorial planes. Second, it ex- 
hibits a hole at its centre. Third, the surface density across 
the disk is non-uniform, with the plasma concentrated into 
two elongated 'clouds' situated at the intersection between 
magnetic and rotational equators. These clouds are best seen 
in the (j) — 180° panel of the figure. In the case of a Ori E, the 
existence of such clouds has been inferred from Ha measure- 



ments (e.g.,|Groote k Hunger] 1982| [Bolton et al.|l987[ ), and 



Townsend et al. (2005) have demonstrated how the clouds 



are simultaneously responsible for the distinctive spectro- 
scopic and photometric variability exhibited by the star. 

To explore further the degree of correspondence with 
the RRM model. Fig. [To] plots the quantities 



and 



R{(Jd) 



i) 



Q"d,RFHD 
C^d,RRM 



(62) 



(63) 



in the x — y magnetic equatorial plane, where the subscripts 
RFHD and RRM denote values obtained from the RFHD sim- 
ulations and the RRM model, respectively. These quantities 



represent the ratio of disk surface densities dd, as defined by 
eqn. (54), and the diff erence in the disk centroid colatitudes 



Od introduced in ^5.2 



The figure reveals that once again the agreement be- 
tween simulations and model is good. The maximum differ- 
ences in surface density are below the 25% level, and - as 
already discussed - are strongly correlated with those disk 
regions that have not undergone any expansion events. The 
differences in centroid colatitude angle do not rise above a 
few degrees; they tend to be most positive around = 0°, 
and most negative around = 180°, so the RFHD disk 
is slightly closer to the magnetic equatorial plane than the 
RRM disk. Based on these findings, we can conclude that the 
disk plasma distribution predicted by the RRM model fur- 
nishes a close approximation to the distribution determined 
via the physically more-sophisticated (yet computationally 
more expensive) RFHD approach. 

Of course, the significant caveat here is that the RRM 
model is limited to consideration only of the cool plasma 
in the disk. Whereas, the RFHD simulations also encom- 
pass the hot post-shock plasma responsible, for instance, 
for magnetospheric X-ray emission. To illustrate the spa- 
tial and thermal distribution of this plasma, we introduce a 
two-temperature emission measure density (EMD), 

^(Ti,T2) = p j n,np5{T-T')dz, 



dT' , (64) 



where 5{) is the Dirac delta function, and Ue is the electron 
number density (cf. eqn. 40). The EMD characterizes the 



radiative recombination emission by plasma in the temper- 
ature range Ti < T < T2. Fig. [11 [ shows images of 8 at the 
end of the simulations, for the 'face on' (0 = 0°, 180°) and 
'edge-on' (0 = 108°) viewing aspects, and for temperature 
ranges that correspond to thermal emission at optical/UV 
(T < 10^ K), extreme UV (10^ K < T < 10^ K), soft X-ray 
(10^ K < T < 10'^ K) and hard X-ray (T > lO'^ K) energies. 

The optical/UV images confirm the concentration of 
Ha-emitting plasma into two clouds. However, at higher 
temperatures a different distribution emerges. In the EUV 
range, the plasma is situated primarily in thin cooling layers 
on either side of the disk. Because the cooling from 10^ K 
to 10^ K is so efficient (on account of the large number of 
metal lines available), these layers are under-resolved in the 
simulations, and appear fragmented into many small islands 
of emission. Similarly fragmented cooling layers are seen in 
the outer parts of the soft X-ray EMD images, but the pre- 
ponderance of the emission in this 10^K<T<10^K range 
comes from a toroidal belt surrounding the star. This belt 
is primarily composed of the hot, dense plasma associated 
with the siphon flows discussed in J 5.2 (and see also Fig.|6|. 

The soft X-ray belt is enclosed by a torus of hard X- 
ray emission, seen in the T > 10^ K images to extend out 
to about 4i^p. Beyond this radius, an asymmetry develops 
on either side of the disk, with those field lines that pass 
near the rotational poles showing much weaker emission 
than those that avoid the poles. This asymmetry is most 
evident in the = 108° (edge-on) panel; observe how the 
emission is strongest in the northern magnetic hemisphere 
on the left-hand side of the image, but vice- versa on the 
right-hand side. The origin of the asymmetry can be seen 
in the = 0° images of Fig. [g] note the significant differ- 
ence in the density of the post-shock regions on either side 
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Figure 9. The proton column density A^p at the end of the simulations, viewed from six different rotational azimuths (f). The oblate star 
is shown in grey at the centre of each panel; arrows indicate the magnetic and rotation axes (the latter being the vertical, fixed one), 
and the curved arcs show field lines having magnetic shell parameters L = 5, 10 and magnetic azimuths ^ = 0°, 60°, 120°, . . . , 300°. 
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Figure 10. The surface density ratio (left) and centroid colatitude difference (right) between the RFHD simulations and the RRM 
model, across the magnetic equatorial (x-y) plane. 
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distribution 

P(T) : 



dlnT 



J J 8(0,T) dxi dyi 



(65) 



Fig. [12] plots P as a function of temperature, for the 
= 0° case shown in Fig. The figure reveals an es- 
sentially bimodal DEM distribution, with a narrow peak at 
T = 22, 500 K = T* corresponding to the cool disk plasma, 
and a broad peak centred at T 4 x 10^ K associated with 
the extended regions of hard X-ray emission seen in the top- 
most images of Fig. |ll| These regions are not subject to any 
appreciable occultation by the star, and as a consequence 
there is almost no change to the X-ray peak seen in Fig. ^] 
as the rotational azimuth 6 is varied. 



Figure 12. The differential emission measure distribution V, 
plotted as a function of temperature T for the cf) = 0° case shown 
in Fig.^^ The ordinate scale is chosen to emphasize the emission 
at EUV and X-ray energies; the truncated peak at low tempera- 
tures extends up to ~ 2 x 10^'-' cm~^. 



of the disk, which translates directly into the 8 asymmetry. 
The density difference itself arises as a result of the position- 
ing of the shock fronts; the disk in the = 0° meridional 
plane is situated in the northern magnetic hemisphere, and 
the shock fronts in this hemisphere are closer to the stellar 
surface, implying a higher density in the post-shock region. 

To complete our discussion of the magnetospheric emis- 
sion, we consider the differential emission measure (DEM) 



6 DISCUSSION 

The analysis in preceding sections is largely directed toward 
comparing the results from the RFHD simulations against 
the corresponding predictions of the RRM formalism. This 
comparison confirms that the two complementary frame- 
works for modeling massive-star magnetospheres are in good 
agreement. However, the RFHD approach has a far broader 
scope than the RRM model, and we have only begun to 
explore its potential applications. Our preliminary investi- 
gation of a cr Ori E-like star has already revealed a variety of 
novel phenomena, such as disk oscillations and centrifugal 
heating. In §6.2| we discuss the prospects for the next logical 
step of comparing RFHD simulations against observations 
of magnetic stars. In §6.3| we review the limitations of the 
RFHD approach, and suggest how these might be overcome 
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in future studies. First, however, we examine the relation- 
ship between RFHD and previous work. 



6.1 Relation to previous work 

The RFHD approach presented in this paper has allowed 
(to our knowledge) the first-ever time-dependent, three- 
dimensional simulation of a massive- star magnet osphere. In 
developing the new approach, we have drawn extensively on 
previous investigations of stellar magnet ospheres. From the 
perspective of disk accumulation, these include the RRM 
model (TO05), and the earlier rigid-field models advanced 
by 'Michel k Sturrock"fT974) and Nakajima { 1985). From the 
perspective of wind dynamics, the key narrative has been the 



MCWS paradigm of Babel & Montmerle (1997a). 



Indeed, the modeling strategy employed by 'B abel <S] 
Montmerle ( |1997a| ) directly foreshadows the RFHD ap- 
proach, and a comparison between the two is appropriate. 
These authors also adopt a rigid-field ansatz, and likewise 
cast the flow in terms of independent 1-D hydrodynamical 
problems that include the effects of rotation and radiative 
cooling. However, they restrict their analysis to the 2-D ax- 
isymmetric case of an aligned dipole (/3 = 0°), and instead 
of constructing global solutions they treat the wind, post- 
shock and disk regions separately, combining them in a post 
hoc step based on matching conditions. Moreover, rather 
than using a hydrodynamical code they treat the flow in the 
wind and post-shock regions by stationary integration of the 
time-independent momentum equation. This rules out any 
possibility of simulating dynamical phenomena such as disk 
oscillations or the fallback of material close to the star. 

For these reasons, we regard the RFHD approach as a 
significant advance beyond the MCWS model that augured 
it. What of other treatments? As discussed in the intro- 
duction, MHD simulations (e.g., ud-Doula Owocki||2002 ) 
tend to be impractical when field lines become almost rigid. 
Therefore, the direct overlap between MHD and RFHD is 
naturally limited - although this may change if the RFHD 
approach can successfully be extended to open field topolo- 
gies, as is discussed below in §6.3.3| Even for stars accessible 
to MHD simulation, the RFHD approach will retain some 
advantage due to its relatively low computational costs. 



6.2 Comparison with observations 

In presenting the results from the 3-D magnetosphere model 
(^, we avoid any specific comparison with observations. 
This reflects the focus of the paper towards introducing 
the RFHD approach and exploring the physical processes at 
work in massive-star magnet ospheres. Here, we outline those 
specific areas where we believe a confrontation between mod- 
els and observations will prove most fruitful. However, apart 
from brief remarks on the correspondence of our results 
to recent and historical observations, we defer the detailed 
quantitative modeling to subsequent papers. 

Townsend et al.| ( [200 5) have already found that the 



RRM model furnishes a good agreement to the Ka and pho- 
tometric variability of the a Ori E. Since the RFHD simula- 
tions predict a disk plasma distribution similar to the RRM 



model (cf. ^5.3), the improvements brought by the simula- 
tions - at optical wavelengths - are likely to be incremental. 



Nevertheless, the radiative transfer treatment adopted by 
[Townsend et al. ( 2005 ) is at the simplest possible level, and 
lacks a realistic treatment of atomic physics. Thus, there is 
ample room for progress in this area. 

The scope for progress is greater at other wavelengths, 
however. In their Chandra survey of M17, [BrQos et al.| ( [2007t 
uncovered a population of B0-B3 stars characterized by hard 
(up to ^ 5keV) X-ray emission. The DEM distributions 
presented in §5.3| are generally consistent with this level of 
hardness, suggesting a magnetic wind-shock origin for the B 
stars' X-rays. To test this hypothesis at a qualitative level, 
the RFHD simulations can be combined with standard emis- 



sion codes such as APEC (Smith et al. 112001), to synthesize 



spectra for direct comparison against the observations. 

Such modeling should also help to understand the X- 
ray emission of known-magnetic Bp stars. Some of these are 
notable on account of the absence of any significant X-ray 



detections (see, e.g. Drake et al. 1994 Czesla & Schmitt 



2007|). In the case of a Ori E X-rays are seen, but they show 



an ambiguous character: a two-temperature fit to Chan- 
dra measurements gives kT ^ 0.7,2.2 — 3.8 keV (Skinner 



et al. 2004), whereas a corresponding fit to XMM-Newton 



data finds kT ^ 0.3, l.lkeV ( Sanz-Forcada, Franciosini & 



[Pallavicini |2004 ) . The RFHD simulations based loosely on 
the parameters of a Ori E predict a non- varying and rather 
harder spectrum, with a DEM distribution peaking around 
kT ^ 3.5 keV (cf. Fig. 12). This harder-than-observed spec- 



trum could be a consequence of our neglect of thermal con- 



duction (cf. ^2.7). Alternatively, it could be that the mass- 



loss rate adopted in the simulations is too high. (A larger 
M leads to a higher-density wind, which in turn pushes 
the shocks further from the star where the flow velocity is 
higher.) In this respect, we note that there is mounting evi- 
dence from a number of independent diagnostics that mass- 
loss rates of OB stars have historically been overestimated 
(see, e.g. Smith Owocki||2006[ and references therein). 

At radio wavelengths, the emission properties of mag- 
netic Bp stars are rather more consistent than for X-rays. 
However, the physical mechanisms responsible remain some- 



what controversial (e.g., Drake 1998). In a recent paper, 



Trigilio et al. ( 2004 ) advance a 3-D model incorporating gy- 



rosynchrotron emission from mildly-relativistic electrons ac- 
celerated in an equatorial current sheet. This current sheet 
lies outside the inner, rigid-field regions of the magneto- 
sphere, and cannot be modeled directly using the RFHD ap- 



proach (although see J 6.3.3). However, Trigilio et al. (2004) 



suggest that the thermal plasma trapped in the post-shock 
cooling regions of the inner magnetosphere plays an impor- 
tant role in modulating (via free-free absorption) the radio 
emission. The RFHD approach will prove useful in testing 
this idea. 

At UV wavelengths, there is good potential for head- 
way to be made. A number of magnetic massive stars were 
observed extensively by lUE (see |Smith Groote|2001 and 
references therein), leading to empirical models for equato- 
rially focused wind streams (e.g.. Shore Brown|1990 ) that 
prefigured the MCWS model. In the case of a Ori E and 



similar He-strong stars, Smith & Groote (2001) find that 



the strengths of UV lines are consistent with absorbing col- 
umn densities of 10^^-10^^ cm . Allowing for the fact that 
these authors assumed a covering factor of 100% (whereas 
magnetospheric disks are in fact closer to 10%), such values 
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are consistent with the A^p 3 x 10^^ cm~^ column densities 
found in the RFHD simulations. The task now is to inves- 
tigate whether the simulations can reproduce the detailed 
absorption profiles measured by lUE. 

On a more speculative note, the RFHD approach may 
be able to shed some light on particle acceleration pro- 
cesses around massive stars. "Bell" (1978a b) and iBlandford 



6.3.2 Rigid field ansatz 

A key ingredient of the RFHD approach is the ansatz that 
field lines remain rigid. This applies as long as the mag- 
netic Lorentz force can balance any competing forces act- 
ing perpendicular to field lines with only minimal distortion 
from the assumed unstressed configuration, taken here to 
be a dipole. In the wind outflow regions, the relevant com- 
petition can be characterized in terms of the global wind 
magnetic confinement parameter, 77*, defined by |ud-Doula| 



& Ostriker ( 1978| ) independently suggested that Fermi ac- 
celeration of particles in astrophysical shocks could ex- 
plain the cosmic-ray energy spectrum up to the 'knee' at 

- 10^ - 10^ GeV. Typically, supernova remnants are con- 

sidered the most likely locations for this to occur fStanev that the Alfven radius - where closed field lines become 

opened into a radial configuration by the wind - scales as 



& Owocki (2002). These authors' MHD simulations show 



2004). However, it seems possible that the circumstellar en- 



vironments of strong-field massive stars could also be sites 
for particle acceleration. Babel k, Montmerle| (p-997a) sug- 
gest that the highly sub-Alfvenic nature of the magneti- 
cally channeled wind shocks (due to the near-rigid field) 
means that the second- order Fermi process should be effec- 
tive. They claim that relativistic electrons produced in this 
manner can explain the radio emission from magnetic mas- 
sive stars, without the need for the current sheets invoked 



by |Trigilio et al.| ( |2004 ). 

We conjecture that the same mechanism operating on 
protons and/or ions (which are far less sensitive to inverse 
Compton losses than electrons) might ultimately lead to the 
production of energetic 7 rays. With recent advent of high- 
sensitivity ground-based Cherenkov telescopes such as HESS 
( Hinton|2004t and VERITAS ( |Holder|2006t , which have the 
sensitivity and angular resolution to detect these 7 rays, the- 
oretical progress on this issue would be particularly timely. 



6.3 Limitations 

To bring the discussion to a close, we briefly review the 
present limitations of the RFHD approach, and (where pos- 
sible) suggest how these limitations might be overcome in 
future extensions to the basic approach. 



6.3.1 Cross-field coupling 



In deriving the adopted expression ( 33 ) for the radiative ac- 
celeration grad, the polar velocity derivative {dv/dO)r is ne- 
glected. Although the approximation is generally quite rea- 
sonable, it necessarily suppresses any coupling between the 
flow on adjacent field lines that may be important in setting 
the scale of knots and cross- field discontinuities (^5.2). 



To include this coupling, we contemplate an extension 
to the basic RFHD approach that involves conducting the 1- 
D simulations lying in the same meridional plane in parallel. 
By this, we mean that each numerical timestep advances the 
flow data of all coplanar field lines together. This way, the 
necessary velocity data to calculate {dv/dr)^ are available 
throughout the simulation, and it it not necessary to make 
the approximation ( 30 ) . The additional computational costs 



of performing the 1-D simulations in parallel are quite rea- 
sonable; beyond the additional memory overhead, the only 
significant issue is that the numerical timestep is limited by 
the Courant criterion as applied to all field lines together. 
This means that the simulations advance at the rate of the 
'slowest' (shortest-Courant time) field line. 



The parameters adopted for the RFHD simulations (cf. 
Table [l]), together with a polar field strength ^ 10^ G 



corresponding to the value inferred for a Ori E (Land- 
street Borra|1978 ) , imply a global confinement parameter 
?7* ?^ 10 . The conclusion is therefore that the wind should 
substantially distort field lines only for radii r > ta ~ 50i^* 
- significantly further out than the maximum shell param- 
eter L — 11.2 considered in our simulations. Although this 
simple analysis does not account for the presence of shocks, 
or for the effects of the Coriolis force, these additional pro- 
cesses can be expected to incur modest (order-unity) modifi- 
cations to the total flow energy. Hence, their inclusion should 
not appreciably modify the basic 77* -based analysis for the 
overall competition between field and flow in regions outside 
the disk. 

However, the radial increase of the centrifugal force 
makes it a potentially important limiting factor in a rigid- 
field approach, particularly for the secularly accumulating 
plasma in the outer regions of the disk. In fact, as discussed 
in the Appendix of TO05, the centrifugal force acting on this 
plasma should eventually become stronger than the avail- 
able magnetic tension, leading to the kind of centrifugally 
driven breakout seen in the MHD simulations of |ud-Doula,| 
Townsend & Owocki ( 2006 ) . The timescale for breakout de- 



creases sharply with local disk radius, asymptotically as r~ ; 
moreover, the overall normalization of this timescale scales 
with a disk confinement parameter that is quite analogous 
to, and has a similar magnitude to, the wind parameter 77* 
(see, e.g., TO05, their eqns. A7 and A8). 

In the RFHD simulations, and again adopting the in- 
ferred field strength for a Ori E, the outermost, L = 11.2 
field line has a disk breakout time of ~ 2 Ms, roughly an 
order of magnitude shorter than the simulation duration. 
Indeed, over this duration all field lines having L > 6 should 
undergo one or more breakout episodes, thus formally vio- 
lating the rigid-field ansatz when applying the RFHD sim- 
ulations of these outer regions to a Ori E. But in practice, 
these regions make only a minor contribution to the disk 



emission (cf. the lower images of Fig. 11 ), implying that the 
model should still be applicable for interpreting optical line 
diagnostics like Ha. Moreover, once field lines reconnect af- 
ter a breakout event, the wind shocks should quickly reform 
to nearly their characteristic asymptotic strength, meaning 
that the associated X-ray emission signatures are also likely 
to be only weakly affected. While further direct MHD sim- 
ulations would be helpful in testing these expectations, it 
seems that models based on the rigid-field ansatz here still 
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provide good approximate predictions of key observational 
diagnostics. 



6.3.3 Open field topologies 

Recall that the general RFHD approach is not necessarily 
limited to any particular magnetic topology, such as the 
dipole form adopted here. It merely requires that this topol- 
ogy be pre-specified and independent of the actual flow. 
For example, the basic approach could even be applied to 
a topology in which a dipole field is opened into a radial 
configuration in the region outside the Alfven radius. Phys- 
ically such opening occurs because of the stresses of wind 
outflow, but phenomeno logically it can be pre-specified by 
invoking current sheets 



(e.g.: 



Connerney et al. 



1981 



Tsy- 



ganenko 1989 Tsyganenko & Peredo 1994) and/or source 
surfaces (e.g., Altschuler &: Newkirk|1969 ) for regions above 
ta- In models of the solar wind, the recent investigation by 



Riley et al. 12006, ) confirms that such source-surface meth- 
ods lead to global field topologies that are in good agreement 
with full MHD solutions. 

In the context of massive stars, corresponding open-field 
RFHD simulations could provide a first attempt toward a 3- 
D model for the effects of the field on open regions of wind 
outflow. Although not contributing appreciably to optical 
or X-ray emission, the inclusion of such outflow regions in 
an RFHD model could prove a good initial basis for syn- 
thesizing phase-dependent UV wind absorption profiles, for 
comparison against archival lUE spectra. 



7 SUMMARY 

We have presented a new Rigid-Field Hydrodynamics ap- 
proach to modeling massive-star magnetospheres. Within 
the rigid-field ansatz, the flow along each field line is treated 
as an independent 1-D hydrodynamical problem. Using the 
VH-1 code, we performed many separate 1-D simulations of 
differing field lines, and pieced them together to build up 
3-D model of a a Ori E-like star. 

These results from these simulations showcase the po- 
tential of the RFHD approach, both by confirming the find- 
ings of previous studies (e..g, wind collision shocks, disc ac- 
cumulation), and by revealing new phenomena (e.g., disk 
oscillations, siphon flows, centrifugal heating). Hence, we an- 
ticipate that the new approach will prove a powerful tool in 
future investigations of massive-star magnetospheres. 
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APPENDIX A: THE MAGNETOHYDRO STATIC 
LIMIT 



magnetohydrostatic equilibrium along the field line. 



dp 
d^ 



d^eff 



(A2) 



Assuming that the stationary plasma has cooled to the stel- 



lar surface temperature T* , the equation of state ( 34 ) is used 
to derive the density distribution 



(A3) 



where C is a constant of integration, and a* is the isothermal 
sound speed introduced in §3.1| 

In the vicinity of a local minimum (as sampled along 
a field line), the effective potential can be approximated by 
the second-order Taylor expansion. 



$eff(s) ^ $0 + 



[s-sof 



(A4) 



where so is the arc distance coordinate of the minimum. 
Here, primes are used to denote differentiation with respect 
to s, while a subscript indicates evaluation at the mini- 
mum s — sq. (The $0 term in the expansion is by definition 
zero.) Substituting this expression into eqn. (A3) reveals a 
Gaussian density distribution. 



with po a constant based on C, and 



ho ; 




(A5) 



(A6) 



the characteristic scale height. 

In the case of a dipole field aligned with the rotation 
axis (i.e., /3 = 0), the effective potential minima are situated 
at the field-line summits, and thus so = 0. It can then be 
shown that 



3GM, 



(A7) 



in the limit r^/R^ (see TO05, their eqn. 18), where 



1/3 



3-Rp 
2w^ 



(A8) 



is the Kepler co-rotation radius where the rotation velocity 
coincides with the orbital velocity. In the same limit, the 
disk scale height is approximated by 



ho ^ ho 



3a* / Rl 



2w \ GM, 



(A9) 



which depends only on the rotation rate and the stellar pa- 
rameters. 



The Rigidly Rotating Magnetosphere (RRM) model is based 
on the magnetohydrostatic limit, where all velocities and 
time derivatives in the Euler equations vanish. In this 

limit, the momentum equation becomes 



dp 
d^ 



■ Pgeff • 



(Al) 



(the grad term is not present because the radiative accel- 
eration scales with the velocity gradient; see J 2.5). With 
the effective gravity expressed in terms of a scalar effective 
potential (cf. eqn. 16), we recover an equation describing 



APPENDIX B: THE RADIAL VELOCITY 
DERIVATIVE 



To evaluate the radial velocity derivative in eqn. (29), we 
first expand it as 



dv 
dL 



dL 
dr 



The right-hand side of this expression includes the velocity 
gradients both along and across the field. To eliminate the 
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cross-field gradient {dv/dL)s, we similarly expand the polar 
velocity derivative as 



ds 



dv 
dL 



dL 



(B2) 



Combining these two expressions, and after some straight- 
forward algebra, we obtain 



ds 



(B3) 



This equation is exact; however, under the assumption (30) 



that the polar velocity derivative can be neglected, the sec- 
ond term on the right-hand side vanishes, so that 



dv\ _ f dv\ 

Or J n 



Here, the second equality follows from eqns. ^ and ([7|, and 
gives the final form (31 ) for the radial velocity derivative. 



m ^ 0. These normal modes have eigenfrequencies set by the 
characteristic equation 



V2m 



ho 



(C5) 



The fundamental (m = 0) mode has a zero frequency, and 
corresponds not to an oscillation, but to the addition of more 
mass to the disk. Higher order modes all conserve disk sur- 
face density [i.e., (crd)p = 0], because 



r 

J —C 



e ^ Hm{x) dx ■■ 



f 

J oo 



e"" Hm{x)Ho{x) dx = (C6) 



for m = 1,2,3,...; the first equality is because Ho{x) = 
1, while the second follows from the orthogonality of the 
Hermite polynomials with respect to t he weighting function 
e~^ (see Abramowitz Stegun|l972 ). 



APPENDIX C: LINEAR PERTURBATION 
ANALYSIS 

To investigate small-amplitude departures from the magne- 
tohydrostatic equilibrium discussed in Appendix |A] we em- 
ploy a linear analysis. Equations ([T]) and ([2| are subjected to 
small-amplitude disturbances; discarding terms of quadratic 
or higher order in the perturbation amplitude leads to the 
system of equations 

^ + i£(^P«p) = 0, (CI) 

Here, the subscript p denotes Eulerian (fixed-position) per- 
turbations, while quantities without subscripts refer to 
the equilibrium state. As before, the grad term has been 
dropped, and the geff term is expressed in terms of the ef- 
fective potential $eff; no term containing the unperturbed 
velocity v appears, because the equilibrium state is taken to 
be at rest. 

To solve these equations, we make the approximations 
that (i) all perturbations share a ti me d ependence of the 



form e^^*; (ii) the Taylor expansion (A4) may be used to 
model the spatial variation of $eff; (hi) the derivative dA/ds 
can be neglected, and (iv) the perturbations are isothermal, 
with the temperature remaining at T* . With these approxi- 
mations, the governing equations for the spatial dependence 
of pp in the vicinity of a potential minimum s = so can be 
reduced to the single second-order equation, 

hlp'; + 2is-so)p'p+(^'^+2^ p^ = 0. (C3) 

As before, primes denote differentiation with respect to s. 



and ho was defined in eqn. (A6). Subject to the boundary 



conditions that pp ^ for (s — so) ±oo, the eigensolutions 
are readily found as 

Pp(s) = am e-(^-^«>'/'^o Hmlis - so)/ho] (C4) 

where 

CLtTi IS a constant setting the amplitude of the pertur- 
bations, and Hm is the Hermite polynomial of integer degree 
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